\chapter{General input module}\label{ch:general}

In this chapter the general structure of the input file for \dalton\
is described, as well as the possible input cards that can be entered in
the \Sec{*DALTON} input module. This input section should always begin an
input file for the \dalton\ program.

\section{General input to DALTON : \Sec{*DALTON}}\label{sec:general}

This input module describes the overall type of calculation that is to
be done. It also contains eight submodules describing
the control of the two different geometry optimization routines,
the control of the three different environment models, PCM, QM3 and PE,
the control of the reading of the molecule and basis set input
  in the ".mol" file,
tuning of the performance of parallel calculations,
as well as control of the general routines for calculating
numerical geometrical derivatives of molecular energies or selected
first- and second-order properties.

We note that \verb|**DALTON| identifies that what follows is general input.
This input section is compulsory to do a calculation, by default \dalton\ does nothing
(will stop with the error message:
{\tt End of file on DALTON.INP, no **DALTON input found}.
Lines before the first line with \verb|**DALTON| are ignored by \dalton .
The general input module stops at the next line starting with two stars: \verb|**anything|.

If more than one node, \dalton\ will automatically perform those modules in parallel
which are implemented in parallel, and Fock matrix calculations will be done AO-direct,
also if you don't specify \verb|.DIRECT|.
In contrast to previous releases of the \dalton\ program, it will not quit in the non-parallelized modules
but run them on the master node while the other nodes are idling.
The keyword \verb|.PARALLEL| of previous releases is still recognized for backward compatibility,
but it is not needed any more.
Note that in order to evaluate the parallelization efficiency\index{parallel efficiency},
a print level of at least 2 is needed in the \Sec{PARALLEL} submodule.

\begin{description}
%\item[\Key{ABACUS}] Alternative keyword to \verb|.PROPERTIES|. Obsolete, do not use.

\item[\Key{CHOLES}] The calculation will use Cholesky decomposed two-electron
integrals~\cite{choint}\index{Cholesky decomposition-based methods} where implemented.
At present, the use of decomposed integrals is only implemented for calculation of energy and
first and second order molecular properties for limited models.
In particular, it is possible to carry out Hartree-Fock
\index{SCF}\index{HF}\index{Hartree--Fock},
Density Functional Theory
\index{DFT}\index{Density Functional Theory},
second-order M{\o}ller-Plesset perturbation theory
(MP2)\index{M{\o}ller-Plesset!second-order}\index{MP2}, and
CC2 coupled cluster calculations\index{CC2}\index{Coupled Cluster}. Note that
the decomposition of differentiated integrals is not implemented yet and, therefore,
gradients cannot be computed. Moreover, Cholesky CCSD(T)~\index{CCSD(T)}
calculations, where the energy dominators are Cholesky decomposed,
are carried out using standard integrals (or in a direct fashion).


\item[\Key{DIRECT}] The calculation is to be done in a direct
manner\index{direct calculation}, that
is, the two-electron integrals\index{two-electron integral!direct} are to be
constructed ``on the fly''
and not written to disc as is the default. This keyword will only work
for SCF\index{SCF}\index{HF}\index{Hartree--Fock} wave functions,
Density Functional Theory
calculation\index{DFT}\index{Density Functional Theory} and for
coupled cluster (CC) calculations\index{CC}\index{Coupled Cluster}. In
HF and DFT calculations the
two-electron integrals (and differentiated
two-electron integrals) will not be written to disc in any part of the
calculation, whereas in the direct CC approach, the two-electron
integrals will be stored in three general-indexed batches, thus
requiring a loop over all basis functions for one of the two-electron
integral indices~\cite{directCC}.

\item[\Key{DOUGLAS-KROLL}] Include scalar relativistic effects by using the
Douglas-Kroll-Hess second-order (DKH2) transformed one-electron potential and kinetic energy Hamiltonian.

\item[\Key{FDE}] activates the use of frozen density embedding (FDE) to include environment effects
(See Chapter~\ref{ch:fde-embedding} for an introduction to FDE). Further options may be specified via
the \Sec{FDE} input section (see subsection~\ref{subsec:fde}).  \index{frozen density embedding} \index{FDE}
\index{environment model} \index{embedding model} \index{multiscale modeling} \index{embedding}
\index{QM/QM} \index{quantum mechanics / quantum mechanics} \index{solvent effects} \index{enviroment effects}
\index{subsystem DFT}


%\item[\Key{HERMIT}] Alternative keyword to \verb|.INTEGRALS|. Obsolete, do not use.

\item[\Key{INPTES}] Test the input of the \Sec{*DALTON INPUT} input
module. The program will abort after the completion of the input test,
and no calculations will be executed.

\item[\Key{INTEGRALS}] Invoke the {\her} and/or the {\eri} module for generating molecular one-
and two-electron integrals.
See Chapter~\ref{ch:hermit} for the \her\ and
Section~\ref{ch:eri}\index{one-electron integral}\index{two-electron integral} for the \eri\ module, respectively.

\item[\Key{ITERATION}]\verb| |\newline
\verb|READ (LUCMD, '(I5)') ITERNR|

Tells the program at which iteration to start the geometry optimization
using the \Sec{*WALK} module\index{iteration number!geometry, start}.
Note that this will {\em not} affect which molecule
input file that is going to be read, as this is handled by the
job script {\tt dalton}. It only determines what number the
output of the predicted molecular geometry will be, as well as where to start
writing information on the files containing information about an IRC
calculation (\verb|DALTON.IRC|) or a dynamical trajectory walk (\verb|DALTON.TRJ|).

\item[\Key{MAX IT}]\verb| |\newline
\verb|READ (LUCMD, '(I5)') ITERMX|

Change the maximum number of geometry iterations\index{iteration number!geometry, max}\index{geometry iterations!max}\index{geometry optimization!max number of iterations}\index{convergence!geometry, max iterations}
that can be done. Default
is~20.  For numerical
differentiation/vibrational averaging, the number of iterations will
be reset to $6N+1$ (where $N$ is the number of nuclei) as the number
of required iterations for these calculations are well
defined. This number has to be increased in Intrinsic Reaction Coordinate
(IRC)\index{IRC}\index{intrinsic reaction coordinate} or dynamical
trajectory studies\index{dynamics}. However, changing this variable will override this reset
option.

%\item[\Key{MINIMIZE}] Alternative keyword to \verb|.OPTIMIZE|. Obsolete,
%do not use.

\item[\Key{FCKTRA}]\verb| |\newline
\index{integral transformation!parallel"} Requests that the parallel
Fock-matrix based two-electron integral transformation routines
by Hans J\o rgen Aagaard Jensen (2016) are used.

\item[\Key{NEWTRA}]\index{integral transformation!"new"} Requests that the
two-electron integral transformation routines of Bj\"{o}rn Roos
should be used during execution of the program, also if less than
256 orbitals. The "new" transformation is always used if more than 255 orbitals.

\item[\Key{NMDDRV}] Calls for a generalized numerical geometry
differentiation. These routines will take advantage of the full
molecular point group in order to minimize the number of points to be
calculated. What order of derivatives and whether any analytical
derivatives are to be used is determined in the \verb|**NMDDRV| input
module (required when \texttt{.NMDDRV} specified).

\item[\Key{NOATMD}] Turns off checking if atoms are too close.

\item[\Key{OPTIMIZE}] Perform a geometry optimization\index{geometry optimization!.OPTIMIZE module}.
The default is a geometry minimization\index{geometry minimization!.OPTIMIZE module}:
an optimization of the molecular geometry to a stationary point with no
negative molecular Hessian eigenvalues
(a local minimum) will be done using the default first-order
methods. However, this may be changed using appropriate
keywords in the submodule \Sec{OPTIMIZE}, and we refer to examples in
the chapter on potential energy surfaces
(Chapter~\ref{ch:geometrywalks}), and subsection~\ref{subsec:minimize}
describing the input cards for the \Sec{OPTIMIZE} submodule for a more
detailed description of possible options.

%% hjaaj Apr 2011: .PARALLEL is not needed any more, Dalton automatically perform the modules
%% in parallel if more than one MPI node is associated.
%\item[\Key{PARALLEL}] Requests that the calculation of two-electron
%integrals\index{two-electron integral}
%is to be done in parallel\index{parallel calculation}. This also
%implies that the calculation is
%done without writing two-electron integrals to disc. This keyword only
%applies to SCF\index{SCF}\index{HF}\index{Hartree--Fock} wave
%functions and DFT calculations\index{DFT}\index{Density Functional Theory}, but all two-electron integral
%evaluations in an SCF calculation will be done parallel as well as the integration of the exchange-correlation functionals and kernels. More details
%about the parallelization strategy in \dalton\ can be found
%in Ref.~\cite{pndjhapdkrthhkcpl253}.

%The keyword requires that the program has been installed and compiled
%with the appropriate preprocessor directives for an MPI
%installation\index{MPI}.
%
%If MPI\index{MPI} is used as message passing\index{message passing}
%interface, no further keywords are
%needed, as the number of nodes will be set equal to the number of
%nodes asked for when submitting the job.
%Note that in order to evaluate the parallelization efficiency\index{parallel efficiency},
%a print level of at least 2 is needed in the \Sec{PARALLEL} submodule.

%% hjaaj Apr 2011: .PARNMD activates PARIO, which is not functional in current version.
%\item[\Key{PARNMD}] Option non-functional. Do not use.

%% hjaaj Apr 2011: Dalton now automatically selects the "new" integral transformation,
%% when more than 255 basis functions.
%\item[\Key{PRESORT}]\index{integral sort} Requests that the
%two-electron integrals should be
%sorted and that the integral transformation routines of Bj\o rn
%Roos should be used during execution of the program.
%The keyword is needed if one attempts an MCSCF, CI or MP2 (run through the \sir\ module and not using the Coupled Cluster module) with more than 255 basis functions.

\item[\Key{PELIB}] Include environment effects by using the polarizable embedding (PE)
or polarizable density embedding models\index{polarizable embedding}\index{PE}\index{environment model}\index{embedding model}\index{embedding}\index{QM/MM}.
This will embed the core molecular system, as defined in the \molinp\ file, in a
PE/PDE potential. See further input options
in the \Sec{PELIB} input section (see subsection~\ref{subsec:peqm}). See
Chapter~\ref{ch:embedding} for an introduction to PE and PDE calculations.

\item[\Key{PRIERR}]\verb| |\newline
\verb|READ (LUCMD, *) IPRERR|
\index{error print level} \index{print level!error}
Reads in the print level that is to be used in the DALTON.ERR
file. Default print level is \verb|IPRUSR+1|.

\item[\Key{PRINT}]\verb| |\newline
\verb|READ (LUCMD, *) IPRUSR|
\index{print level!general}
Reads in the print level that is to be used the rest of the subsequent
calculations. Default is a print level of 0.

\item[\Key{PROPERTIES}] Invoke the \aba\ module for the evaluation of static
and dynamic properties, using a previously converged and saved wave function and saved integral files. See Chapter~\ref{ch:abacus}.

\item[\Key{QFIT}] Invoke the calculation of atomic partial charges fitted to reproduce
the molecular electrostatic potential after the wave function has converged.

\item[\Key{RESPONSE}]
Invoke the \resp\ module for the evaluation of static and dynamic
properties, using a previously converged and saved wave function and saved integral files. See Chapter~\ref{ch:response}.

\item[\Key{RUN ALL}]
Invoke all the modules {\her}, {\sir}, {\resp}, and {\aba} for a single point calculation.

\item[\Key{RUN PROPERTIES}]
Invoke the modules {\her}, {\sir}, and {\aba} for a single point calculation.

\item[\Key{RUN RESPONSE}]
Invoke all the modules {\her}, {\sir}, and {\resp} for a single point calculation.

\item[\Key{RUN WAVE FUNCTIONS}]
Invoke the modules {\her} and {\sir} for a single point energy calculation.

%\item[\Key{RUNABA}] Alternative keyword to \verb|.RUN PROPERTIES|. Obsolete, do not use.

\item[\Key{RUNERI}] Force the use of the vectorized 2-electron integral code {\eri} where possible.

%\item[\Key{RUNRSP}] Alternative keyword to \verb|.RUN RESPONSE|. Obsolete, do not use.

%\item[\Key{RUNSIR}] Alternative keyword to \verb|.RUN WAVE FUNCTIONS|. Obsolete, do not use.

\item[\Key{TOTSYM}] Consider only totally symmetric
perturbations\index{symmetry!only totally symmetric perturbations}.
This option only affects geometric perturbations calculated using the
second-order based \Key{WALK} option and static
electric-field perturbations requested through the keyword \Key{POLARI}.

\item[\Key{VECLEN}]\verb| |
\newline
\verb|READ(LUCMD,*) IVECLN|

Set the number of Fock matrices to be used during Fock-matrix constructions in
direct calculations. This function is only of interest for vector machines.
The default is 128. The larger the number, the more memory will be required
in the calculation.

\item[\Key{WALK}]
Do a geometry walk\index{geometry walk}\index{geometry optimization!.WALK module}.
If no input is given in the
\Sec{WALK} input submodule, an optimization of the molecular
geometry to a stationary point with no negative molecular Hessian eigenvalues (a
local minimum) will be done using a second-order method with
analytical molecular Hessians. However, this may be changed by appropriate
keywords in the submodule \Sec{WALK}, and we refer to examples in
the chapter on potential energy surfaces
(Chapter~\ref{ch:geometrywalks}), and subsection~\ref{sec:abawalk}
describing the input cards for the \Sec{WALK} submodule for a more
detailed description of possible options.

\item[\Key{WAVE FUNCTIONS}]
\index{SCF}\index{HF}\index{Hartree--Fock}\index{MP2}\index{MCSCF}\index{CC}
\index{NEVPT2}\index{GASCI}\index{DFT}
Invoke the {\sir} module for the evaluation of SCF, MP2, Coupled
Cluster, MCSCF, NEVPT2, GASCI wave functions as well as DFT calculations, according to input.
Necessary one-electron integrals \em{must} already be evaluated and available (AOPROPER and AOONEINT files).
The two-electron integrals will be generated if needed and not available (AOTWOINT file).
Invoke the modules {\her} and {\sir} for a single point energy calculation using old saved integral files (AOONEINT, AOPROPER, and - if not .DIRECT or .CHOLESKY - AOTWOINT).
See Chapter~\ref{chap:sirius-inpref}.
\end{description}

%hjaaj Feb 2003: next is redundant, do not mention but keep in code for now
%\subsection{End of General input: \Sec{END OF}}
%
%The last input card for the general input section (\Sec{*DALTON})
%may be \Sec{END OF}.

\subsection{Geometry optimization module 1: \Sec{OPTIMIZE}}\label{subsec:minimize}

This submodule is the usual driver for geometry optimizations\index{geometry optimization!.OPTIMIZE module},
but we note that the \Sec{WALK} module contains another second-order
geometry optimization driver for HF, DFT, and MCSCF.
Only for analytical second-order optimization of HF, DFT, and MCSCF can you choose
any of the two geometry optimization modules.
For the unique capabilities of the other module, see subsection~\ref{sec:abawalk}.

The \Sec{OPTIMIZE} module contains both first and second-order methods
for locating minima and transition states (geometry optimization).
\index{first-order geometry optimization}\index{geometry optimization!first-order}
\index{second-order geometry optimization}\index{geometry optimization!second-order}
Most of the molecular Hessian updating\index{molecular Hessian!Hessian update} schemes were taken from
ref.~\cite{thkrprt95} and \cite{Fletcher}.
The implementation of redundant internal coordinates
\index{redundant internal coordinates}\index{coordinate system!redundant internal coordinates}
follows the work of Peng {\it et al.\/}~\cite{cppyahbsmjfjcc17}.
In addition to this, several keywords for VRML visualization are included~\cite{VRML}.

For wave function models where analytical molecular gradients and/or Hessians are not implemented,
the displacements for evaluation of numerical gradients and Hessians
are controlled by \Key{DISPLA} in module \Sec{*NMDDRV} (see subsection\ref{sec:nmddrv}).

\begin{description}

\item[\Key{1STORD}]
Use default first-order method\index{first-order geometry optimization}\index{geometry optimization!first-order}.
This means that the BFGS update will\index{BFGS update}\index{Molecular Hessian!Hessian update!BFGS}
be used, and that the optimization is carried out in redundant internal
coordinates\index{redundant internal coordinates}. Same effect as the
combination of the two keywords \Key{BFGS} and \Key{REDINT}. Since the
\Key{BFGS} method ensures a positive definite Hessian, the
\Key{BOFILL} optimization method is used by default in case of
searches for transition states.

\item[\Key{2NDORD}]
Use default second-order method\index{second-order geometry optimization}\index{geometry optimization!second-order}.
Molecular Hessians will be calculated at every
geometry. The level-shifted Newton method and
Cartesian coordinates\index{Cartesian coordinates} are used. Identical
to specifying the keywords \Key{NEWTON} and \Key{CARTES}.

\item[\Key{BAKER}]
Activates the convergence criteria of Baker \cite{Baker}. The minimum
is then said to be found when the largest element of the gradient
vector (in Cartesian or redundant internal coordinates) falls below
$3.0\cdot 10^{-4}$ and either the energy change from the last
iteration is less than $1.0\cdot 10^{-6}$ or the largest element of
the predicted step vector is less $3.0\cdot 10^{-4}$.

\item[\Key{BFGS}]
Specifies the use of a first-order method\index{first-order geometry optimization}
with the Broyden-Fletcher-Goldfarb-Shanno (BFGS)
update\index{BFGS update}\index{Hessian update!BFGS}
formula for optimization. This is the
preferred first-order method for minimizations, as this update is able
to maintain a positive definite Hessian. Note that this also makes it
unsuitable for transitions state optimization (where one negative
eigenvalue is sought).

\item[\Key{BFGSR1}]
Use a linear combination of the BFGS and the symmetric rank one
updating schemes in the same fashion as Bofill's update. Only suitable
for minimizations.

\item[\Key{BOFILL}]
Bofill's update\cite{jmbjcc15}\index{Bofill's update}\index{Molecular Hessian!Hessian update!Bofill's update}
is the default updating scheme for
transition state optimizations. It's a linear combination of the
symmetric rank one and the PSB updating schemes, automatically giving
more weight to PSB whenever the rank one potentially is numerically
unstable.

\item[\Key{CARTES}]
Indicates that Cartesian coordinates\index{Cartesian coordinates}\index{coordinate system!Cartesian coordinates}
should be used in the optimization. This is the default for
second-order methods.

\item[\Key{CMBMOD}]
Uses a combination of the BFGS update and the model molecular Hessian (diagonal
in redundant internal coordinates). The two have equal weight in the
first iteration of the geometry optimization, then for each subsequent
iteration the weight of the model Hessian is halved. Only suitable for
minimizations.

\item[\Key{CONDIT}]\verb| | \newline
\verb|READ (LUCMD,*) ICONDI|

Set the number of convergence criteria\index{convergence criteria!geometry}
that should be fulfilled before
convergence occurs. There are three different convergence thresholds,
one for the energy\index{geometry convergence criteria!energy change}, one for the gradient
norm\index{geometry convergence criteria!norm of gradient} and one for the step
norm\index{geometry convergence criteria!norm of step}.
The possible values for this variable is therefore between 1 and
3. Default is 2. The three convergence thresholds can be adjusted with
the keywords \Key{ENERGY}, \Key{GRADIE} and \Key{STEP T}.

\item[\Key{CONSTRAINT}]\verb| |\newline
\verb|READ (LUCMD, *) NCON|\newline
\verb|DO I = 1, NCON|\newline
\verb|   READ(LUCMD,*) ICON|\newline
\verb|   ICNSTR(ICON) = 1|\newline
\verb|END DO|

Request a constrained geometry optimization. Only works when using
redundant internal coordinates. The number of primitive coordinates
that should be frozen has to be specified (\verb|NCON|), then a list
follows with the individual coordinate numbers. The coordinate numbers
can be found by first running \dalton\ with the \Key{FINDRE}
keyword. Any number of bonds, angles and dihedral angles may be
frozen. NOTE: Symmetry takes precedence over constraints, if you
{\it e.g.\/} want to freeze just one of several symmetric bonds, symmetry
must be lowered or switched off.

\item[\Key{DELINT}]
Use delocalized internal coordinates\index{delocalized internal coordinates}\index{coordinate system!delocalized internal coordinates}.
These are built up as
non-redundant linear combinations of the redundant internal
coordinates. Performance is more or less the same as for the redundant
internals, but the transformation of displacements (step) is slightly
less stable.

\item[\Key{DFP}]
Specifies that a first-order\index{first-order geometry optimization} method with the
Davidon-Fletcher-Powell (DFP) update\index{DFP update}\index{molecular Hessian!Hessian update!DFP}
formula should be used for optimization. May be used for both
minimizations and transition state optimizations.

\item[\Key{ENERGY}]\verb| | \newline
\verb|READ(LUCMD,*) THRERG|

Set the geometry convergence threshold for the energy (in a.u.). This is one of the three
convergence thresholds\index{convergence!geometry, criteria} (the keywords
\Key{GRADIE} and \Key{STEP T}
control the other two). Default value is the maximum of $1.0\cdot
10^{-6}$ and two times the threshold for the wave function gradient.

\item[\Key{FINDRE}]
Determines the redundant internal coordinate system then quits without
doing an actual calculation. Useful for setting up constrained
geometry optimizations, where the numbers of individual primitive
internal coordinates are needed.

\item[\Key{FREEZE}] \verb| |\newline
\verb|READ (LUCMD, *) NFREEZ|\newline
\verb|DO I = 1, NFREEZ|\newline
\verb|   READ(LUCMD,*) IFREEZ(I)|\newline
\verb|END DO|

Request the \verb|NFREEZ| specified atoms to be frozen in the geometry optimization.
This option enforces cartesian geometry optimization
and disables any use of translation and rotation invariance.
This option was implemented to freeze capping atoms in geometry optimization of
the QM part of a QM/PE system, using polarizable embedding for the environment.
It can of course also be used if the user wants to freeze (clamp) atoms for other reasons.
See also keyword \Key{FRZITR}. Default: all atoms included in geometry optimization.

\item[\Key{FRZITR}]\verb| | \newline
\verb|READ(LUCMD,*) ITRFRZ|

Stop freezing positions of atoms frozen with \Key{FREEZE} from geometry iteration no. \verb|ITRFRZ|.

\item[\Key{GDIIS}]
Use the Geometrical DIIS\cite{pcppjms114}\index{geometrical DIIS}
algorithm to control the
step. Works in much the same way as DIIS for wave functions. However,
the rational function and level-shifted Newton methods are generally
more robust and more efficient. Can only be used for minimizations.

\item[\Key{GEOANA}]
Enables an analysis of the molecular geometry in terms of bond lengths
and bond angles at each new geometry predicted during the optimization
procedure.

\item[\Key{GRADIE}]\verb| |
\newline
\verb|READ(LUCMD,*) THRGRD|

Set the convergence threshold for the norm of the molecular gradient (in
a.u.). This is one of
the three convergence thresholds (the keywords \Key{ENERGY} and
\Key{STEP T} control the other two). Default value is the maximum of
$1.0\cdot 10^{-4}$ and two times the threshold for the wave function
gradient.

\item[\Key{GRDINI}]
Specifies that the molecular Hessian\index{molecular Hessian!reinitialization} should be
reinitialized every time the norm
of the gradient is larger than norm of the gradient two iterations
earlier. This keyword should only be used when it's difficult to
obtain a good approximation to the molecular Hessian during optimization. Only
applies to first-order methods\index{first-order geometry optimization}.

\item[\Key{HELLMA}]
Use gradients and Hessians calculated using the Hellmann-Feynman
approximation. Currently not working properly

\item[\Key{HESFIL}]
Specifies that the initial molecular Hessian
\index{initial Hessian!first-order geometry optimization}\index{molecular Hessian!initial}
should be read from the file \verb|DALTON.HES|. This applies to first-order
methods, and the Hessian in the file must have the correct
dimensions. This option overrides other options for the initial Hessian.

Each time a Hessian is calculated or updated\index{molecular Hessian!DALTON.HES file},
it's written to this file (in Cartesian coordinates). If an
optimization is interrupted, it can be restarted with the
last geometry and the molecular Hessian in \verb|DALTON.HES|, minimizing the
loss of information. Another useful possibility is to transfer
the molecular Hessian from a calculation on the same molecule with another
(smaller) basis and/or a cheaper wave function. Finally, one can go in
and edit the file directly to set up a specific force field.

\item[\Key{INIMOD}]
\index{model Hessian}\index{molecular Hessian!initial!model Hessian}
Use a simple model Hessian~\cite{rlabgkpamcpl241} diagonal in redundant
internal coordinates as the initial Hessian. All diagonal elements are
determined based on an extremely simplified molecular mechanics model,
yet this model provides Hessians that are good starting points for
most systems, thus avoiding any calculation of the exact Hessian. This
is the default for first-order methods.

\item[\Key{INIRED}]
Specifies that the initial molecular Hessian\index{molecular Hessian!initial}
should be diagonal in redundant internal coordinates. The different diagonal
elements are set equal to 0.5 for bonds, 0.2 for angles and 0.1 for
dihedral angles, unless \Key{INITEV} has been specified. If the
optimization is run in Cartesian coordinates, the diagonal internal
molecular Hessian is transformed to Cartesians. Only applies to first-order
methods\index{first-order geometry optimization}.

\item[\Key{INITEV}]\verb| | \newline
\verb|READ(LUCMD,*) EVLINI| \index{molecular Hessian!initial!diagonal}

The default initial molecular Hessian for first-order
minimizations\index{first-order geometry optimization} is the
identity matrix when Cartesian coordinates are used, and a diagonal
matrix when redundant internal coordinates are used. If \Key{INITEV}
is used, all the diagonal elements (and therefore the eigenvalues) are
set equal to the value EVLINI. This option only has effect when
first-order methods are used and \Key{INITHE} and \Key{HESFIL} are
non-present.

\item[\Key{INITHE}]
\index{molecular Hessian!initial}
Specifies that the initial molecular Hessian should be
calculated (analytical molecular Hessian), thus yielding a first step that is
identical to that of second-order methods. This provides an excellent starting
point for first-order methods, but should only be used when the
molecular Hessian can be calculated within a reasonable amount of time. It has only
effect for first-order methods and overrides the keywords
\Key{INITEV} and \Key{INIRED}. It has no effect when \Key{HESFIL} has
been specified.

%\item[\Key{INTERA}]
%Specifies a interactive run. The energy, gradient
%norm\index{norm of gradient}, step length\index{norm of step} and
%molecular Hessian index\index{molecular Hessian!index} is then written to the unit LUERR(=0)
%in each iteration, allowing easy monitoring of the optimization.

\item[\Key{LINE S}]
Turns on line searching\index{line search!geometry optimization}\index{geometry optimization!line search},
using a quartic polynomial. By default this
is turned off, as there seems to be no gain in efficiency. Can only be
used for minimizations.

\item[\Key{M-BFGS}]
A list of old geometries and gradients are kept. At each new point,
displacements and gradient difference for the last few steps are
calculated, and all of these are then used to sequentially update the
molecular Hessian, the most weight being given to the last displacement and
gradient difference. Each update is done using the BFGS formula, and
it's thus only suitable for minimizations. Only applies to first-order
methods\index{first-order geometry optimization}.

\item[\Key{M-PSB}]
This identical to \Key{M-BFGS}, except the PSB formula is used for the
updating. Only applies to first-order methods\index{first-order
geometry optimization}, but it can be used for both minimizations and saddle
point optimizations.

\item[\Key{MAX IT}]\verb| | \newline
\verb|READ(LUCMD,*) ITRMAX|

Read the maximum number of geometry iterations\index{geometry iteration}\index{geometry optimization!iterations}.
Default value is 25.

\item[\Key{MAX RE}]\verb| | \newline
\verb|READ(LUCMD,*) MAXREJ|

Read maximum number of rejected steps\index{rejected geometry step}\index{geometry optimization!rejected step}
in each iterations, default is 3.

\item[\Key{MODE}]\verb| | \newline
\verb|READ(LUCMD,*) NSPMOD|

Only has effect when doing saddle point optimizations\index{mode following}\index{geometry optimization!mode following}.
Determines which molecular Hessian eigenmode should be maximized (inverted in the image
method). By default this is the mode corresponding to the lowest
eigenvalue, {\it i.e.\/} mode 1. If an optimization does not end up at
the correct transition state, it may be worthwhile following other modes
(only the lower ones are usually interesting).

\item[\Key{MODHES}]
Determine a new model molecular Hessian (see \Key{INIMOD})
at every geometry without doing any updating. The model is thus used
in much the same manner as an exact molecular Hessian, though it is obviously
only a relatively crude approximation to the analytical molecular Hessian.

\item[\Key{NEWTON}]
Specifies that a second-order Newton method should be used for\index{geometry optimization!second-order}
optimization---that is, the analytical molecular Hessian will be calculated at
every geometry. By default the level-shifted trust region method will
be used, but it is possible to override this by using one of the two
keywords \Key{RF} or \Key{GDIIS}.

\item[\Key{NOAUX}]
Only has effect when using redundant internal coordinates. The default
for minimizations is to add auxiliary bonds between atoms that are up
to two and half times further apart then regular (chemical)
bonds. This increases the redundancy of the coordinate system, but
usually speeds up the geometry optimization slightly. \Key{NOAUX}
turns this off. For saddle point optimizations and constrained
geometry optimization this is off by default (cannot be switched on).

\item[\Key{NOBREA}]
Disables breaking of symmetry\index{symmetry!breaking}\index{geometry optimization!symmetry breaking}.
The geometry will be optimized within the given symmetry, even if a non-zero molecular Hessian
index is found. The default is to let the symmetry be broken until
a minimum is found with a molecular Hessian index\index{molecular Hessian!index} of
zero. This option only has effect when second-order methods are used.

\item[\Key{NODIHE}]
Only has effect when using redundant internal coordinates. Removes all
coordinates that are dihedral angles, leaving only bonds and
angles. Not too useful, but may be used if one wants to limit the
number of internal coordinates. Constrained geometry optimizations can
sometimes benefit from having all dihedral angles removed (assuming no
dihedral angles needs to be frozen).

%\item[\Key{NOREIN}]
%When first-order methods are used, the
%molecular Hessian\index{molecular Hessian!reinitialization} is reinitialized every
%time the molecular Hessian index\index{molecular Hessian!index} becomes non-zero (due to negative
%eigenvalues). This guarantees that the molecular Hessian describes a minimum,
%but valuable information gathered in the molecular Hessian may be
%lost.
%
%\Key{NOREIN}
%Disables this reinitialization, relying on the
%optimization method to restore the molecular Hessian to its correct structure
%(by locating the area near the minimum). This option is particularly
%useful in conjunction with the keyword \Key{INITHE}, as it is usually
%not meaningful to calculate the molecular Hessian at the initial geometry, then
%resetting it to the identity matrix because some negative eigenvalue
%showed up.

\item[\Key{NOTRUS}]
Turns off the trust radius, so that a full Newton step
is taken in each iteration. This should be used with caution, as
global convergence is no longer guaranteed. If long steps are desired,
it is safer to adjust the initial trust radius and the limits for the
actual/predicted energy ratio.

\item[\Key{PREOPT}]\verb| |\newline
\verb|READ (LUCMD,*) NUMPRE|\newline
\verb|DO I = 1, NUMPRE|\newline
\verb|   READ (LUCMD,*) PREBTX(I)|\newline
\verb|END DO|

First we read the number of basis sets\index{basis sets!preoptimization} that should be
used for preoptimization\index{preoptimization!basis sets}, then we read those
basis set names as strings. These
sets will be used for optimization in the order they appear in the
input. One should therefore place the smaller basis at the
top. After the preoptimization, optimization is performed with the
basis specified in the molecule input file.

\item[\Key{PRINT}]\verb| |
\newline
\verb|READ (LUCMD,*) IPRINT|

Set print level for this module.  Read one more line containing print
level. Default value is 0, any value higher than 12 gives debugging
level output.

\item[\Key{PSB}]
Specifies that a first-order method\index{first-order geometry optimization} with the
Powell-Symmetric-Broyden (PSB)\index{PSB update}\index{molecular Hessian!update!PSB}
update formula should be used for optimization.

%\item[\Key{QUADSD}]
%Not implemented.

\item[\Key{RANKON}]
Specifies that a first-order method with the rank one update\index{rank one update}\index{molecular Hessian update!rank one}\index{molecular Hessian update!MS}\index{molecular Hessian update!SR1}
formula should be used for optimization. This updating is also
referred to as symmetric rank one (SR1) or Murtagh-Sargent (MS).

\item[\Key{REDINT}]
Specifies that redundant internal coordinates\index{redundant internal coordinates}\index{coordinate system!redundant internal coordinates}
should be used in the optimization. This is the default for
first-order methods\index{first-order optimization}.

\item[\Key{REJINI}]
Specifies that the molecular Hessian should be
reinitialized\index{molecular Hessian!reinitialization} after every
rejected step\index{rejected geometry step}, as a rejected step
indicates that the molecular Hessian models the
true potential surface poorly. Only applies to first-order
methods\index{first-order geometry optimization}.

\item[\Key{REMOVE}]\verb| |\newline
\verb|READ (LUCMD, *) NREM|\newline
\verb|DO I = 1, NREM|\newline
\verb|   READ(LUCMD,*) IREM|\newline
\verb|   ICNSTR(IREM) = 2|\newline
\verb|END DO|

Only has effect when using redundant internal coordinates.
Specifies internal coordinates that should be removed. The input is
identical to the one for \Key{CONSTRAINT}, that is one has to specify
the number of coordinates that should be removed, then the number of
each of those internal coordinates. The coordinate numbers can first
be determined by running with \Key{FINDRE} set.

Removing certain coordinates can sometimes be useful in speeding up
constrained geometry optimization, as certain coordinates sometimes
``struggle'' against the constraints. See also \Key{NODIHE}.

\item[\Key{RF}]
Use the rational function method~\cite{abnajsrsjpc89} instead of
\index{rational function}
level-shifted Newton which is the default. The RF method is often
slightly faster than the level-shifted Newton, but also slightly less
robust.

For saddle point optimizations there's a special partitioned rational
function method (used automatically when both \Key{RF}
and \Key{SADDLE} are set). However, this method is both slower and
less stable than the default trust-region level-shifted image method
(which is the default).

\item[\Key{SADDLE}]
Indicates that a saddle point optimization should be performed rather\index{geometry optimization!transition state}\index{transition state}
than a minimization. The default method is to calculate the molecular Hessian
analytically at the initial geometry, then update it using Bofill's\index{Bofill's update}\index{molecular Hessian update!Bofill's update}
update. The optimization is performed in redundant internal
coordinates and using the trust-region level-shifted image method to
control the step. That is by default all the keywords \Key{INITHE},
\Key{BOFILL} and \Key{REDINT} are already set, but this can of course
be overridden by specifying other keywords. If locating the desired
transition state is difficult, and provided analytical molecular Hessians are
available, it may sometimes be necessary to use the \Key{NEWTON}
keyword so that molecular Hessians are calculated at every geometry.

\item[\Key{SCHLEG}]
Specifies that a first-order method\index{first-order geometry optimization}
with Schlegel's updating scheme\index{Schlegel update!geometry optimization}
\cite{Schlegel} should be used. This makes use of all previous
displacements and gradients, not just the last, to update the
molecular Hessian.

\item[\Key{SP BAS}]\verb| | \newline
\verb|READ(LUCMD,*) SPBSTX|

Read a string containing the name of a basis set. When the geometry
has converged, a single-point energy will be calculated using this
basis set\index{basis set!converged geometry}.

\item[\Key{STABILIZE}]\verb| | \newline
\verb|READ(LUCMD,*) ISTBLZ|

Tries to ``stabilize'' the predicted new molecular geometries (and
thus reduce the risk of symmetry breaking) by ignoring all numbers
appearing in the Cartesian coordinates of the atoms beyond digit
number \verb|ISTBLZ|.

\item[\Key{STEEPD}]
Specifies that the first-order steepest descent\index{first-order geometry optimization}
\index{steepest descent} method should be used.
No update is done on the molecular Hessian, so the optimization will be
guided by the gradient alone. The ``pure'' steepest descent method is
obtained when the molecular Hessian is set equal to the identity matrix. Each
step will then be the negative of the gradient vector, and the
convergence towards the minimum will be extremely slow. However, this
option can be combined with other initial molecular Hessians in Cartesian or
redundant internal coordinates\index{Cartesian coordinates}\index{redundant internal coordinates},
giving a method where the main feature is the lack of molecular Hessian updates (static molecular Hessian).

\item[\Key{STEP T}]\verb| | \newline
\verb|READ(LUCMD,*) THRSTP|

Set the geometry convergence\index{convergence!geometry, criteria} threshold for
the norm of the geometry step (in a.u.).
This is one of the three convergence thresholds (the keywords \Key{ENERGY} and
\Key{GRADIE} control the other two). Default value is $1.0\cdot
10^{-4}$.

\item[\Key{SYMTHR}]\verb| | \newline
\verb|READ(LUCMD,*) THRSYM|

Determines the gradient threshold (in a.u.) for breaking of the
symmetry\index{symmetry!breaking}. That is, if the index of the
molecular molecular Hessian\index{molecular Hessian!index of molecular} is non-zero when the gradient norm
drops below this value, the symmetry is broken to avoid unnecessary
iterations within the wrong symmetry. This option only applies to
second-order\index{second-order geometry optimization} methods and when the
keyword \Key{NOBREA} is not present. The default value of this
threshold is $5.0\cdot 10^{-3}$.

\item[\Key{TR FAC}]\verb| | \newline
\verb|READ(LUCMD,*) TRSTIN, TRSTDE|

Read two factors that will be used for increasing and decreasing the
trust radius\index{trust radius!geometry optimization}, respectively, in geometry optimization. Default values are
1.2 and 0.7.

\item[\Key{TR LIM}]\verb| | \newline
\verb|READ(LUCMD,*) RTENBD, RTENGD, RTRJMN, RTRJMX|

Read four limits for the ratio between the actual and predicted
energies. This ratio indicates how good the step is---that is, how
accurately the quadratic model describes the true energy
surface. If the ratio is below \verb|RTRJMN| or above
\verb|RTRJMX|, the step is rejected. With a ratio between
\verb|RTRJMN| and \verb|RTENBD|, the step is considered bad an the
trust radius\index{trust radius!geometry optimization} decreased to less than the step
length. Ratios between \verb|RTENBD| and \verb|RTENGD| are
considered satisfactory, the trust radius is set equal to the norm
of the step. Finally ratios above \verb|RTENGD| (but below
\verb|RTRJMX|) indicate a good step, and the trust radius is given
a value larger than the step length. The amount the trust radius
is increased or decreased can be adjusted with \Key{TR FAC}. The
default values of \verb|RTENBD|, \verb|RTENGD|, \verb|RTRJMN| and
\verb|RTRJMX| are 0.4, 0.8, -0.1 and 3.0 respectively.

\item[\Key{TRSTRG}]
Specifies that the level-shifted trust region method should be used to
control the step. This is the default, so the keyword is actually
redundant at the moment. Alternative step control methods are \Key{RF}
and \Key{GDIIS}.

\item[\Key{TRUSTR}]\verb| | \newline
\verb|READ(LUCMD,*) TRSTRA|

Set initial trust radius\index{trust radius!geometry optimization} for the geometry optimization. This
will also be the
maximum step length for the first iteration. The trust radius is
updated after each iteration depending on the ratio between predicted
and actual energy change. The default trust radius is 0.5 a.u.

\item[\Key{VISUAL}]
Specifies that the molecule should be
visualized\index{visualization}\index{VRML}\index{geometry!visualization with VRML},
writing a VRML file of the molecular geometry. {\em{No optimization
will be performed when this keyword is given}}. See also related
keywords \Key{VR-BON}, \Key{VR-COR}, \Key{VR-EIG} and \Key{VR-VIB}.

\item[\Key{VRML}]
Specifies that the molecule should be
visualized\index{visualization}\index{VRML}. VRML files describing
both the initial and final geometry
will be written (as \verb|initial.wrl| and \verb|final.wrl|). The file
\verb|final.wrl| is  updated in each iteration, so that it always
reflects the latest geometry. See also related keywords \Key{VR-BON},
\Key{VR-COR}, \Key{VR-EIG} and \Key{VR-VIB}.

\item[\Key{VR-BON}]
Only has effect together with \Key{VRML} or \Key{VISUAL}. Specifies
that the VRML files should include bonds\index{bonded atoms} between nearby
atoms\index{visualization}\index{VRML}. The
bonds are drawn as grey cylinders, making it easier to see the
structure of the molecule. If \Key{VR-BON} is omitted, only the
spheres representing the different atoms will be drawn.

\item[\Key{VR-COR}]
Draws $x$-, $y$- and $z$-axis in the VRML scenes with
geometries. Somewhat useful if one is struggling to build a reasonable
geometry by adjusting coordinates manually.

\item[\Key{VR-EIG}]
Only has effect together with \Key{VRML} or
\Key{VISUAL}\index{visualization!eigenvectors}\index{VRML}.
Specifies that the eigenvectors of the molecule (that is the
eigenvectors of the molecular Hessian, which differs from the normal modes as
they are not mass-scaled) should be visualized. These are written to
the files \verb|eigv_###.wrl|.

\item[\Key{VR-SYM}]
Draws in all symmetry elements of the molecule as vectors (rotational
axes) and semi-transparent planes (mirror planes).

\item[\Key{VR-VIB}]
Similar to \Key{VR-EIG}, but more useful as it draws the actual normal
mode vectors (the mass-weighted eigenvectors). These are written to
the files \verb|norm_###.wrl|. Keyword only has effect when a
vibrational analysis has been requested.

\end{description}


\subsection{Parallel calculations : \Sec{PARALLEL}}

This submodule controls the performance of the parallel
version\index{parallel calculation}
of \dalton . The implementation has been described
elsewhere~\cite{pndjhapdkrthhkcpl253}. \dalton\ only supports MPI as
message passing interface in the current release.

\begin{description}

\item[\Key{DEBUG}]\verb| |
Transfers the print level from the master\index{slave!debug}
to the slaves, otherwise the print level on the slaves will always be
zero. Only for debugging purposes.

\item[\Key{DEGREE}]\verb| |\newline
\verb|READ (LUCMD,*) NDEGDI|

Determines the percent of available tasks\index{parallel tasks} that
is to be distributed in
a given distribution of tasks, where a distribution of tasks is defined
as the process of giving batches to {\em all} slaves. The default is
5\% , which ensures that each slave will receive 20 tasks during one
integral evaluation, which will give a reasonable control with the
idle time of each slave.

\item[\Key{PRINT}]\verb| |\newline
\verb|READ (LUCMD, *) IPRPAR|

Specify the print level for the parallel control routines in a parallel calculation.
A print level of at least 2 is needed in order to be able to evaluate the
parallel efficiency\index{parallel efficiency}. A complete
timing for all nodes will be given if the print level is 4 or higher.
\end{description}

\subsection{Polarizable embedding model: \Sec{PELIB}}
\label{subsec:peqm} \index{polarizable embedding}\index{PE}\index{environment model}\index{embedding model}\index{embedding}\index{QM/MM}\index{multiscale model}\index{multiscale}\index{quantum mechanics / molecular mechanics}\index{solvent effects}\index{enviroment effects}

This input section controls calculations using the polarizable
embedding (PE)~\cite{pemodel1,pemodel2}
and polarizable density embedding (PDE)~\cite{pde1,pde2,pde3} models
which can include effects from a structured
environment on a core molecular system described using quantum mechanics.
See Chapter~\ref{ch:embedding} for an introduction to
calculations using the PE and PDE models.

\begin{description}

\item[\Key{POTENTIAL}]\verb| |\newline
\verb|READ (LUCMD, *) POTFILE|\verb| |\newline
This option can be used to specify a non-default name of the \potinp\ input file that contains the embedding potential parameters. The default name is \verb|POTENTIAL.INP|.

\item[\Key{PDE}]\verb| |\newline
\verb|READ (LUCMD, *) HDF5FILE|\verb| |\newline
This option can be used to specify a non-default name of the \verb|hdf5| input file that contains the PDE operators. The default name is \verb|standard.h5|.

\item[\Key{DIRECT}]\verb| |\newline
Use the direct solver to determine the induced dipole moments. It will explicitly build a classical response matrix of size $3S(3S+1)/2$, where $S$ is the number of polarizable sites and is therefore only recommendable for smaller molecular systems. Note that this solver is not parallelized. The default is to use the iterative solver.

\item[\Key{ITERATIVE}]\verb| |\newline
\verb|READ (LUCMD, *) THRITER| (optional)\verb| |\newline
Use the iterative solver to determine the induced dipole moments. This is the default. The convergence threshold defaults to $1.0\cdot10^{-8} > |\mu^{[k]} - \mu^{[]k-1]}|$, where $\mu$ is a vector containing all induced dipoles and $k$ is the iteration index, but can also be provided with this option.

\item[\Key{BORDER}]\verb| |\newline
1) \verb|READ (LUCMD, *) BORDER_TYPE, RMIN, AUORAA|\verb| |\newline
2) \verb|READ (LUCMD, *) BORDER_TYPE, REDIST_ORDER, RMIN, AUORAA, NREDIST|\verb| |\newline
Controls the handling of the border between the core molecular system
and its environment described by the embedding potential. There are
two mutually exclusive schemes: 1) \verb|BORDER_TYPE = REMOVE| and 2)
\verb|BORDER_TYPE = REDIST|. The first option will remove all
multipoles and polarizabilities that are within the given distance
\verb|RMIN| from any atom in the core molecular system. The
\verb|AUORAA| variable specifies whether the distance threshold is
given in \angstrom{} (AA) or \bohr{} (AU). The second option will
redistribute parameters that are within the given threshold
\verb|RMIN| from any atom in the core system to nearest sites in the
environment. The order of multipoles up to which will be redistributed
is determined by the \verb|REDIST_ORDER| variable,
e.g. \verb|REDIST_ORDER = 1| means that only charges will be
redistributed and all other parameters removed,
\verb|REDIST_ORDER = 2| means charges and dipoles are redistributed
and so on. The sign of \verb|REDIST_ORDER| specifies if the
polarizabilities are redistributed. Positive means that the
polarizabilities are removed and negative means redistributed. The
number of sites that parameters on a given site are redistributed to
is determined by the \verb|NREDIST| variable which can be between 1
and 3. The default is to redistribute charges within 2.2 \bohr{} to its
nearest site and removing all other parameters.

\item[\Key{DAMP}] (deprecated)\verb| |\newline
\verb|READ (LUCMD, *) DAMP| (optional)\verb| |\newline
\textbf{This option has been deprecated and will be removed in a future release (see below for other options).} Damp the interactions between induced dipole moments using Thole's exponential damping scheme~\cite{thole_damp_1,thole_damp_2}. The default damping coefficient is the standard 2.1304~\cite{thole_damp_2}.

\item[\Key{DAMP INDUCED}]\verb| |\newline
\verb|READ (LUCMD, *) IND_DAMP| (optional)\verb| |\newline
Damp the electric fields from induced dipole moments using Thole's exponential damping scheme~\cite{thole_damp_1,thole_damp_2}. The default damping coefficient is the standard 2.1304~\cite{thole_damp_2}.

\item[\Key{DAMP MULTIPOLE}]\verb| |\newline
\verb|READ (LUCMD, *) MUL_DAMP| (optional)\verb| |\newline
Damp the electric fields from permanent multipole moments using Thole's exponential damping scheme~\cite{thole_damp_1,thole_damp_2}. This option requires polarizabilities on all sites with permanent multipole moments. The default damping coefficient is the standard 2.1304~\cite{thole_damp_2}.

\item[\Key{DAMP CORE}]\verb| |\newline
\verb|READ (LUCMD, *) CORE_DAMP| (optional)\verb| |\newline
\verb|READ (LUCMD, *) NALPHAS|\verb| |\newline
\verb|DO I = 1, NALPHAS|\verb| |\newline
\verb|    READ (LUCMD, *) ISO_ALPHA(I)|\verb| |\newline
\verb|END DO|\verb| |\newline
Damp the electric fields from the electrons and nuclei in the core region based on Thole's exponential damping scheme~\cite{thole_damp_1,thole_damp_2}. The damping coefficient is optional unless isotropic polarizabilities are present. The default damping coefficient is the standard 2.1304~\cite{thole_damp_2}. Standard polarizabilities from ref~\citenum{thole_damp_2} have been implemented and will be used if none are given as input. However, only H, C, N, O, F, S, Cl, Br and I are available.

\item[\Key{GSPOL}]\verb| |\newline
Activate the ground-state polarization approximation, i.e.\ freeze the embedding potential according to the ground-state density. This means that the polarizable environment is self-consistently relaxed during the optimization of the ground-state density/wave function of the core molecular system and then kept frozen in any following response calculations.

\item[\Key{NOMB}]\verb| |\newline
Remove many-body effects in the environment. This is done by deactivating interactions between induced dipole moments.

\item[\Key{EEF}]\verb| |\newline
Request the calculation of effective dipole integrals, i.e.,
include the effective external field effect in dipole property calculations.
NOTE: Activating this option will affect all dipole properties, including the dipole moment
of the core quantum region.

\item[\Key{PBC}]\verb| |\newline
Use periodic boundary conditions (PBCs) in the calculation of induced dipoles. This
option requires that the simulation box is defined in the input potential file (see Chapter~\ref{ch:embedding} for details).

\item[\Key{RESTART}]\verb| |\newline
Use any existing files to restart calculation.

\item[\Key{CUBE}]\verb| |\newline
1) \verb|READ (LUCMD, *) STD_GRID|\verb| |\newline
2) \verb|READ (LUCMD, *) OPTION|\verb| |\newline
2) \verb|  IF (OPTION == 'GRID') THEN|\verb| |\newline
2) \verb|     READ (LUCMD, *) XSIZE, XGRID, YSIZE, YGRID, ZSIZE, ZGRID|\verb| |\newline
2) \verb|  END IF| \verb| |\newline
\verb|  IF (OPTION == 'FIELD') THEN|\verb| |\newline
\verb|      FIELD = .TRUE.| \verb| |\newline
\verb|  END IF| \verb| |\newline
Create cube file for the core molecular system containing the
electrostatic potential due to the final converged polarizable
embedding potential. The grid density can be specified either
using 1) standard grids (\verb|COARSE| (3 points/\bohr{}), \verb|MEDIUM| (6 points/\bohr{})
or \verb|FINE| (12 points/\bohr{}) and in all cases 8.0 \bohr{} are added in each direction) or 2) the \verb|GRID|
option which gives full control of the cube size and density, and
requires an additional input line specifying the extent added (in \bohr{}) in each direction and
density (in points/\bohr{}). The default grid density is \verb|MEDIUM| and default cube size
is the extent of the molecule plus 8.0 \bohr{} in plus and minus each
Cartesian coordinate.
If the \verb|FIELD| option is given then also three cube files
containing the Cartesian components of the electric field from the
embedding potential will be created.

\item[\Key{FMM}]\verb| |\newline
Activate FMM for evaluations of static and induced multipole fields.

\item[\Key{THETA}]\verb| |\newline
\verb|READ (LUCMD, *) THETA|\verb| |\newline
Value of the multipole acceptance criterion for FMM.

\item[\Key{EXPANSION ORDER}]\verb| |\newline
\verb|READ (LUCMD, *) EXPANSION_ORDER|\verb| |\newline
Multipole expansion order used in FMM.

\item[\Key{NCRIT}]\verb| |\newline
\verb|READ (LUCMD, *) NCRIT|\verb| |\newline
Maximum number of particles in an FMM box (before it is split).

\item[\Key{INT\_SCHEME}]\verb| |\newline
\verb|READ (LUCMD, *) INT_SCHEME|\verb| |\newline
Approximate coupling scheme to use (if any). Choices are EXACT (no approximation, default), SINGLE\_CENTER (single-center multipole expansion), ESPF (multicenter ESP-fitted multipole expansion), and FMM\_FOCK (FMM-like environment-centered expansion).

\item[\Key{INT\_ORDER}]\verb| |\newline
\verb|READ (LUCMD, *) INT_ORDER|\verb| |\newline
Order of the multipole expansion (ESPF or SINGLE\_CENTER) or environment local expansion (FMM\_FOCK) used in the approximate environment coupling.

\item[\Key{INT\_THETA}]\verb| |\newline
\verb|READ (LUCMD, *) INT_THETA|\verb| |\newline
Relevant only for FMM\_FOCK. Parameter controlling near-field/far-field splitting for the FMM\_FOCK coupling method. Typical values are somewhere between 0.0 and 1.0. Lower values the parameter places more interactions in the near-field.

\item[\Key{INT\_RCUT}]\verb| |\newline
\verb|READ (LUCMD, *) INT_RCUT|\verb| |\newline
Relevant only for SINGLE\_CENTER and ESPF. Cutoff radius (in bohrs) beyond which the approximate coupling is used. Larger values of the radius place more interactions in the near-field.

\item[\Key{ESPF\_GRID\_TYPE}]\verb| |\newline
\verb|READ (LUCMD, *) ESPF_GRID_TYPE|\verb| |\newline
Relevant only for ESPF. Type of grid to use when building the ESP-fitting equations. Options are SOLVENT or VDW. 
The SOLVENT grid uses the positions of the environment sites as grid points (as specified in the potential file), provided they are within a certain cutoff radius of the quantum region INT\_RFIT.
The VDW grid uses a molecular VDW surface as grid points. The grid is created using the method by Saff and Kuijlaars.

\item[\Key{INT\_RFIT}]\verb| |\newline
\verb|READ (LUCMD, *) INT_RFIT|\verb| |\newline
Relevant only for ESPF with the SOLVENT ESPF\_GRID\_TYPE. Radius (in bohrs) within which solvent sites are used in the ESP-fitting. 

\item[\Key{ESPF\_SOLVER}]\verb| |\newline
\verb|READ (LUCMD, *) ESPF_SOLVER|\verb| |\newline
Relevant only for ESPF. Type of numerical solver to use in ESPF when solving the ESP-fitting equations. Choices are SPTRS (linear equation solver) or SVD (singular value decomposition).

\item[\Key{SVDTOL}]\verb| |\newline
\verb|READ (LUCMD, *) SVDTOL|\verb| |\newline
Relevant only for ESPF with the SVD numerical solver. Tolerance below which singular values in the SVD are discarded.

\item[\Key{NF\_INCORE}]\verb| |\newline
\verb|READ (LUCMD, *) NF_INCORE|\verb| |\newline
Keep near-field interaction integrals in memory (default for approximate couplings).

\item[\Key{NO\_NF\_INCORE}]\verb| |\newline
\verb|READ (LUCMD, *) NO_NF_INCORE|\verb| |\newline
Do not keep near-field interaction integrals in memory.

\item[\Key{SKIPMUL}]\verb| |\newline
Skip calculation of multipole--multipole interaction energy (default)

\item[\Key{NOSKIPMUL}]\verb| |\newline
Do not skip calculation of multipole--multipole interaction energy.

\item[\Key{ISOPOL}]\verb| |\newline
Converts all anisotropic polarizabilities into isotropic polarizabilities.

\item[\Key{ZEROPOL}]\verb| |\newline
Remove all polarizabilities.

\item[\Key{ZEROMUL}]\verb| |\newline
\verb|READ (LUCMD, *) ZEROMUL_ORDER| (optional) \verb| |\newline
Remove all multipoles of order \verb|ZEROMUL_ORDER| and up. The default is to remove all dipoles and higher-order multipoles (i.e. \verb|ZEROMUL_ORDER| = 1).

\item[\Key{VERBOSE}]\verb| |\newline
Verbose output. Currently this will print the final converged induced dipole moments.

\item[\Key{DEBUG}]\verb| |\newline
Debug output. Prints the total electric field and the induced dipole moments in each iteration. WARNING: for large systems this will produce very large output files.

\item[\Key{OLDFIELD}]\verb| |\newline
Activates the legacy implementation of direct multipole field evaluation. For debugging purposes.

\end{description}

\subsection{QM/MM model: \Sec{QM3}}
\label{sec:qm3}

Not documented. Examples can be found by grepping for "QM3" in the "DALTON/test" directory in the distribution. See also \Key{PELIB} and \Sec{PELIB}
(introduction in Chapter~\ref{ch:embedding}) for a newer more general and efficient implementation with many of the same capabilities.

\subsection{Polarizable continuum model: \Sec{PCM}}
\label{subsec:pcm}

This input section, together with the \Sec{PCMCAV} input section,
controls the polarizable continuum model (PCM) to describe the
environment effects. See Section~\ref{sec:pcm} for an introduction to
calculations using PCM in {\dalton}\index{PCM}.

\begin{description}

\item[\Key{SOLVNT}]\verb| |\newline
\verb|READ (LUCMD, *) RZSOL|\verb| |\newline
Solvent. The following solvents are supported in
{\dalton}

\begin{tabular}[h]{l|l}
Formula  & Keyword \\
\hline
H2O      & WATER \\
CH3OH    & METHANOL \\
C2H5OH   & ETHANOL \\
CHCL3    & CLFORM \\
CH2CL2   & METHYLCL \\
C2H4CL2  & 12DICLET \\
CCL4     & TETRACLC \\
C6H6     & BENZENE \\
C6H5CH3  & TOLUENE \\
C6H5CL   & CLBENZ \\
CH3NO2   & NITROMET \\
C7H16    & N-EPTANE \\
C6H12    & CYCLOHEX \\
C6H5NH2  & ANILINE \\
CH3COCH3 & ACETONE \\
THF      & TETHYDFU \\
DMSO     & DIMETSOX \\
CH3CN    & ACETONIT \\
\hline
\end{tabular}\newline

It is also possible to define a solvent by the keywords \Key{EPS},
\Key{EPSINF} and \Key{RSOLV}.

\item[\Key{NEQRSP}]\verb| |\newline
Non-equilibrium contributions of the solvent to the response
calculation.


\item[\Key{ICESPH}]\verb| |\newline
\verb|READ (LUCMD, *) ICESPH|\verb| |\newline
Define how to put the spheres when generating the PCM cavity. ICESPH
can either be 0, 1 or 2.
\begin{itemize}
\item[0] automatic assignment of spheres on all atoms following the
  tabulated radii values. This is the default.
\item[1] reading the sphere centers coordinates in the \molinp\ file,
  as (X,Y,Z,R), with R the radius of the sphere. Be aware that this
  option is not working if molecular symmetry is detected.
\item[2] specifying the nuclei associated with the sphere. This is
  used together with the \Key{NESFP} keyword and, from the
  \Sec{PCMCAV} section (see below), the \Key{INA} and \Key{RIN}
  keywords
\end{itemize}

\item[\Key{NESFP}]\verb| |\newline
\verb|READ (LUCMD, *) NESFP|\verb| |\newline
Number of spheres (only if ICESPH is defined as 2).

\item[\Key{EPS}]\verb| |\newline
\verb|READ (LUCMD, *) EPS|\verb| |\newline
Defines the static dielectric constant.

\item[\Key{EPSINF}]\verb| |\newline
\verb|READ (LUCMD, *) EPSINF|\verb| |\newline
Defines the optical dielectric constant.

\item[\Key{RSOLV}]\verb| |\newline
\verb|READ (LUCMD, *) RSOLVI|\verb| |\newline
Defines the solvent radius.

\item[\Key{PRINT}]\verb| |\newline
\verb|READ (LUCMD, *) IPRPCM|\verb| |\newline
Print level for the PCM calculation (default 0).

\end{description}

\subsection{The PCM cavity: \Sec{PCMCAV}}
\label{subsec:pcmcav}

This input section is only used right after the *PCM section
(Section~\ref{subsec:pcmcav}), and is used to specify the cavity.

\begin{description}

\item[\Key{AREATS}]\verb| |\newline
\verb|READ (LUCMD, *) AREATS|\verb| |\newline
Defines the maximum area of a tessera in the initial tesselation of a
sphere. (default 0.4)

\item[\Key{INA}]\verb| |\newline
\verb|DO I = 1, NESFP|\newline
\verb|  READ (LUCMD, *) INA(I)|\verb| |\newline
\verb|END DO|\newline
List of atoms to put a sphere on

\item[\Key{RIN}]\verb| |\newline
\verb|DO I = 1, NESFP|\newline
\verb|  READ (LUCMD, *) RIN(I)|\verb| |\newline
\verb|END DO|\newline
The radii of each sphere, given in \angstrom{}.

\item[\Key{ALPHA}]\verb| |\newline
\verb|DO I = 1, NESFP|\newline
\verb|  READ (LUCMD, *) ALPHA(I)|\verb| |\newline
\verb|END DO|\newline
Scaling factor for the sphere radii (default 1.2). If parameter
ALPHA(I) is greater than 0 the I-th radius RIN(I) is multiplied by
ALPHA(I). This allows to use radii=R(van der Waals)*ALPHA in the
calculation of electrostatic (and eventually SCF disp-rep)
contribution, and radii = R(van der Waals) for the cavitation.

\item[\Key{CENTER}]\verb| |\newline
\verb|DO I = 1, NESFP|\newline
\verb|  READ (LUCMD, *) XE(I),YE(I),ZE(I)|\verb| |\newline
\verb|END DO|\newline
xyz coordinates of the center of cavity spheres, given in \angstrom{}.

\end{description}


\subsection{Frozen density embedding model: \Sec{FDE}}
\label{subsec:fde}
\index{frozen density embedding}
\index{FDE}
\index{environment model}
\index{embedding model}
\index{multiscale modeling}
\index{embedding}
\index{QM/QM}
\index{quantum mechanics / quantum mechanics}
\index{solvent effects}
\index{enviroment effects}
\index{subsystem DFT}

This input section controls calculations using the frozen density embedding model.
See Chapter~\ref{ch:fde-embedding} for more details on the method and implementation details.

\begin{description}

\item[\Key{EMBPOT}]\verb| |\newline
\verb|READ (LUCMD, '(A60)') name |\verb| |\newline
Specifies the name of the file containing the embedding potential over an integration grid (used to construct
the matrix representation of the embedding potential), which is subsequently added to the one-electron Hamiltonian.
This option can be used to specify a non-default name of the file that contains the embedding potential. The default name is \verb|EMBPOT|.

%\item[\Key{GRIDOU}]\verb| |\newline
%\verb|READ (LUCMD, *) name|\verb| |\newline
%Enables the calculation and export of the total electrostatic potential and electron density, gradient and Hessian to the file
%with name \verb|name|. The outputted file will be in XML format.
%
%\item[\Key{LEVEL}]\verb| |\newline
%\verb|READ (LUCMD, *) ex_level|\verb| |\newline
%Defines the wavefunction with which to calculate the total electrostatic potential and electron density, gradient and Hessian
%which will be exported. Default is \verb|HF|, supported options are: \verb|HF|, \verb|DFT|, \verb|MP22|, \verb|CC2|, \verb|CCSD|, \verb|MCSCF|.
%
%\item[\Key{EXONLY}]\verb| |\newline
%\verb|READ (LUCMD, *) name|\verb| |\newline
%Enables the calculation and export of the total electrostatic potential and electron density, gradient and Hessian to the file
%with name \verb|name|, while disabling the import of embedding potentials. The outputted file will be in XML format.
%
\item[\Key{PRINT}]\verb| |\newline
\verb|READ (LUCMD, *) pl|\verb| |\newline
Sets the print level for the FDE contributions. The default is \verb|0|.

\end{description}


\subsection{Potential derived charges: \Sec{QFIT}}
\label{subsec:qfit}

Directives controlling the evaluation of charges (and higher ranks of multipole moments, see \verb|MPRANK| below) fitted to reproduce the molecular electrostatic potential evaluated on a grid surrounding the molecule. Currently, the only built-in grid is a Connolly-Merz grid where \verb|NSHELL| shells around the molecule are constructed from van der Waals spheres placed on each nuclei of the molecule. The initial spheres are scaled by \verb|VDWSCL| and succesive spheres are scaled in increments of \verb|VDWINC|. Alternatively, a file with grid points given through the \verb|MEPFIL| keyword can be provided instead.

\begin{description}
\item[\Key{MPRANK}]\verb| |\newline
\verb|READ (LUCMD,*) MPRANK|

Specifies the highest rank of the multipole moments which are fitted to reproduce the electrostatic potential. It is possible to specify either 0 (charges), 1 (charges and dipoles) or 2 (charges, dipoles and quadrupoles). Please note, that using the dipole constraint (see \verb|CONSTR| below) with a rank $>=1$ is meaningless and will be ignored. The default value is 0 (charges).

\item[\Key{CONSTR}]\verb| |\newline
\verb|READ (LUCMD,*) CONSTR|

Set the type of constraint on the derived charges. Options are \verb|NONE|, \verb|CHARGE| or \verb|DIPOLE| which sets no constraints, constrains the sum of potential derived charges to be equal to the molecular charge or constrains the sum and dipole moment of the potential derived charges to be equal to the molecular charge and permanent molecular dipole moment, respectively. Setting the constraint to \verb|NONE| will lead to additional charge being generated or removed and should be used carefully. The default value is \verb|CHARGE|.

\item[\Key{VDWSCL}]\verb| |\newline
\verb|READ (LUCMD,*) VDWSCL|

Sets the scaling factor, $R_{scl}$, for the van der Waals radii of the nuclei when generating the surface. The default value is 1.4.

\item[\Key{VDWINC}]\verb| |\newline
\verb|READ (LUCMD,*) VDWINC|

Sets the increment of scaling, $R_{inc}$, for successive shells past the first. See option \verb|NSHELL|. The total scaling factor for each van der Waals sphere for shell $n$ can be calculated from $R = R_{vdw}\cdot (R_{scl}+(n-1)\cdot R_{inc})$. The default value is 0.2 \angstrom{} per shell.

\item[\Key{NSHELL}]\verb| |\newline
\verb|READ (LUCMD,*) NSHELL|

Generate \verb|NSHELL| shells of points around the molecule. The default value is 4.

\item[\Key{PTDENS}]\verb| |\newline
\verb|READ (LUCMD,*) PTDENS|

The surface point density. Default value is 0.28 points au$^{-2}$.

\item[\Key{MEPFIL}]\verb| |\newline
\verb|READ (LUCMD,*) MEPFIL|

Specifies the file to load from the work directory containing coordinates of a molecular surface.
The format for the surface is an \verb|.xyz|-file with a title specifying either \angstrom{} (AA) or \bohr{} (AU) units:
\begin{verbatim}
20
AA
X    -2.916   0.154  -0.000
X    -3.683   0.221  -0.000
.      .       .       .
.      .       .       .
.      .       .       .
X     1.291   2.156  -0.000
X     1.687   2.819   0.001
\end{verbatim}

\end{description}

%Not documented yet. Examples can be found by grepping for "PCM" in the "DALTON/test" directory in the distribution.

\subsection{Geometry optimization module 2: \Sec{WALK}}
\label{sec:abawalk}

Directives controlling one of the two second-order
geometry\index{second-order geometry optimization}
optimizations  as well as the
execution of dynamical walks and numerical
differentiation\index{numerical differentiation} in
calculations of Raman intensities and optical activity\index{Raman
  intensity}\index{Raman optical activity}\index{ROA},
appear in the \Sec{WALK} section.

\begin{description}
\item[\Key{ANHARM}]
Requests that a determination of the cubic force field is to be determined.
By default this is done calculating numerical derivatives of analytical
molecular Hessians in Cartesian coordinates.

\item[\Key{DISPLA}]\verb| |\newline
\verb|READ (LUCMD, *) DISPLC|

Displacement taken in a numerical differentiation\index{numerical
differentiation!geometry}. This applies both for a numerical molecular
molecular Hessian\index{molecular Hessian!numerical}, as well as in calculation of Raman
intensities and optical activity\index{Raman intensity}\index{Raman optical activity}\index{ROA}.
Read one more line specifying value.  Default is~10$^{-4}$~bohr.
However,
note that this variable does not determine the displacements used
when evaluating numerical gradient for use in first-order geometry
optimizations controlled by the \Sec{OPTIMIZE} module, for example with
MP2 or CI wave\index{MP2}\index{CI}\index{M{\o}ller-Plesset!second-order}\index{Configuration Interaction} functions,
these displacements are controlled by the \Key{DISPLA} keyword in the \Sec{*NMDDRV} module.

\item[\Key{DYNAMI}]
Perform a ``dynamic walk''\index{dynamics}: integrate the
classical equations of motion\index{equation of motion} for the nuclei
analytically on a locally
quadratic surface. The method is discussed in
Ref.~\cite{theuhjajcpl173} as well as in Section~\ref{sec:dynamic}.

\item[\Key{ECKART}]\verb| |\newline
\verb| DO I = 1, NUCDEP|\newline
\verb|    READ (LUCMD,'(7X,F17.10,2F24.10)')|\newline
\verb|&        (ECKGEO(J,I), J = 1, 3)|\newline
\verb| END DO|

During a vibrational averaging, ensure that the properties are
transformed to the appropriate Eckart axis system.  The coordinate
system given should be that of the equilibrium geometry of the
molecule. The coordinates should be given in \bohr{}, and should be given
for all symmetry generated atoms in the order given in the input.

\item[\Key{EIGEN}] Take a step to the boundary of the trust
region\index{trust radius}
along the eigenvector mode\index{eigenvector} specified by \Key{MODE}.

\item[\Key{FRAGME}]\verb| |\newline
\verb|READ (LUCMD, *) NIP|\newline
\verb|READ (LUCMD, *) (IPART(IP), IP = 1, NIP)|

Identify which fragments\index{molecular fragments} atoms belong to in
a dynamic walk.  Read one more line specifying the number of
atoms (the total number of atoms in the molecule), then one more
line identifying which fragment an atom belongs to. The atoms in the
molecule are given a number, different for each fragment. See also the
discussion in Sec.~\ref{sec:dynamic}.

\item[\Key{GRDEXT}] Perform a gradient extremal-based\index{gradient extremal} optimization. The algorithm used in this kind of optimization is
thoroughly described in Ref.\cite{pjhjajthtca73}. This is the default walk
type if the  index of the critical point searched is higher than
1. See also the discussion in Sec.~\ref{sec:gradext}.

\item[\Key{HARMON}]\verb| |\newline
\verb|READ (LUCMD, *) ANHFAC|

Threshold for harmonic dominance.  Read one
more line specifying value. Default is~100. This is  another
way of changing the criterion for changes of the trust
radius\index{trust radius}. See also the keyword \Key{TRUST}.

\item[\Key{IMAGE}] Locate a transition state\index{transition state} using a
trust-region-based image surface\index{image surface} minimization.
Note that only a
point with a molecular Hessian index of~1 can currently be located with this
method, not higher-order stationary points. See also the discussion in
Sec.~\ref{sec:image}.

\item[\Key{INDEX}]\verb| |\newline
\verb|READ (LUCMD,*) IWKIND|

Desired molecular Hessian index\index{molecular Hessian!index} (strictly speaking, of the
totally symmetric block of the molecular Hessian) at the optimized geometry.
Read one more line specifying value.  Default is~0 (minimum).
Note that a stationary point with the wrong molecular Hessian index will not
be accepted as an optimized geometry.

\item[\Key{IRC}]\verb| |\newline
\verb|READ (LUCMD, *) IRCSGN|

Set the geometry walk to be  an
Intrinsic Reaction Coordinate (IRC)\index{intrinsic reaction coordinate}\index{IRC} as described in
Ref.~\cite{kfacr14}. Read one
more line containing the sign (-1 or 1) of the reaction coordinate. It
cannot be decided in advance which reaction pathway a specific sign is
associated with. See also the discussion in Sec.~\ref{sec:irc}.

%\item[\Key{ISOTOP}]\verb| |\newline
%\verb|READ (LUCMD, *) NIS|\newline
%\verb|READ (LUCMD, *) (ISOTPS(IS), IS = 1, NIS)|
%
%Specify the isotopic constitution\index{isotopic constitution} of the
%molecule under investigation.
%This is most interesting in dynamic walks\index{dynamics},
%when using mass-scaled atomic coordinates\index{mass-weighted coordinates}.
%Note that the internal structure of
%\aba\ uses Cartesian coordinates, and for vibrational analysis alone
%the keyword \Key{ISOTOP} in the \Sec{VIBANA} input section is to be
%used. Note also that for defining center of mass
%coordinates or for vibrational averaging, there is a similar \verb|.ISOTOP| keyword in the general
%input module, and this will become the only keyword for specifying
%isotopic substitutions in later versions of the program. The
%\Key{ISOTOP} keyword in this input module (\Sec{WALK}) is to become
%obsolete.

\item[\Key{KEEPSY}] Ensure that the symmetry of the molecule is not
broken. The threshold for determining a mode as breaking symmetry is
controlled by the keyword \Key{ZERGRD}.

\item[\Key{MASSES}] Mass-scale the atomic
coordinates\index{mass-weighted coordinates}.  This is the
default for dynamic walks\index{dynamics}, gradient
extremal\index{gradient extremal} walks and in calculations
of Intrinsic Reaction Coordinates
(IRCs)\index{intrinsic reaction coordinate}\index{IRC}.

\item[\Key{MAXNUC}]\verb| |\newline
\verb|READ (LUCMD, *) XMXNUC|

Maximum displacement\index{displacement of atom} allowed for any one
atom as a result of the geometry update.  Read one more line
specifying value.  Default is~0.5.

\item[\Key{MAXTRU}]\verb| |\newline
\verb|READ (LUCMD, *) TRUMX1|

Set the maximum arc length in an Intrinsic Reaction Coordinate (IRC)
walk\index{intrinsic reaction coordinate}\index{IRC}. Read one more
line containing the maximum arc length.
Default is~0.10. Note that this arc length is also affected by the
\Key{TRUST} keyword, and if both are specified, the arc length will
be set to the minimum value of these to.

\item[\Key{MODE}]\verb| |\newline
\verb|READ (LUCMD,*) IMODE|

Mode to follow in level-shifted Newton optimizations for transition
states\index{transition state}.  Read one more line specifying
mode. Default is to follow
the lowest mode (mode~1).

\item[\Key{MODFOL}] Perform a mode-following (level-shifted
Newton) optimization. This is the default for minimizations and
localization of transition states. See also discussion in
Section~\ref{sec:modfol}.

\item[\Key{MOMENT}]\verb| |\newline
\verb|   READ (LUCMD, *) NSTMOM|\newline
\verb|   DO IP = 1, NSTMOM|\newline
\verb|      READ (LUCMD, *) ISTMOM(IP), STRMOM(IP)|\newline
\verb|   END DO|

Initial momentum for a dynamic walk\index{dynamics}\index{momentum}.
Read one more line specifying
the number of modes to which there is added an initial momentum. Then
read one line for each of these modes, containing first the number of
the mode, and then the momentum. The default is to have no momentum.
See also the section describing how to perform a dynamic
walk, Sec.~\ref{sec:dynamic}.

\item[\Key{NATCON}] Use the natural connection\index{natural connection} when orthogonalizing
the predicted molecular orbitals at the new geometry. By default the
symmetric connection is used.

\item[\Key{NEWTON}] Use a strict
Newton--Raphson\index{Newton-Raphson step} step to update
the geometry. This means that no trust region will be used.

\item[\Key{NO CENTRIFUGAL FORCES}] Do not include contributions from
centrifugal forces when calculating vibrationally averaged geometries
at a finite temperature.

\item[\Key{NOGRAD}]\verb| |\newline
\verb|READ (LUCMD, *) NZEROG|\newline
\verb|READ (LUCMD, *) (IZEROG(I), I = 1, NZEROG)|

 Set some gradient elements  to zero.  Read
one more line specifying how many elements to zero, then one
or more lines listing their sequence numbers.

\item[\Key{NOORTH}] The predicted molecular orbitals at the new
geometry are {\em not} orthogonalized. Default is that the orbitals
are orthogonalized with the symmetric connection. Orthogonalization
can also be done with the natural
connection~\cite{joklbkrthpjtca90}. See the keyword \Key{NATCON}.

\item[\Key{NOPRED}] No prediction of the energy of the wave function
at the updated geometry.

\item[\Key{NORMAL}] Do the calculation of effective (vibrationally
averaged) geometries in normal coordinates. This will restrict the
calculation of the effective geometry to one isotopic species (by
default the most abundant one).

\item[\Key{NUMERI}] Do a numerical
differentiation\index{numerical differentiation}, for instance when
calculating Raman intensities or Raman optical activity, see
Sections~\ref{sec:ramanint} and~\ref{sec:vroa}.

\item[\Key{PRINT}]\verb| |\newline
\verb|READ (LUCMD,*) IPRWLK|

Set the print level in the prediction of new geometry steps.  Read one
more line containing print level. Default value is the value of
\verb|IPRDEF| in the general input module.

\item[\Key{RATLIM}]\verb| |\newline
\verb|READ (LUCMD, *) RTRMIN, RTRGOD, REJMIN, REJMAX|

Limits on ratios between predicted and
observed energy change.  Read one more line specifying four
values~(*).  These are respectively the bad prediction ratio, good
prediction ratio, low rejection ratio and high rejection ratio.
Defaults are~0.4, 0.8, 0.1, and 1.9.

\item[\Key{REJECT}] Signals that the previous geometry step was
rejected\index{rejected geometry step}, and the trust
region\index{trust radius} is
reduced. This keyword is
used in case of restarts to tell the program that when the program
was stopped, the last geometry was in fact rejected.

\item[\Key{REPS}]\verb| |\newline
\verb|READ (LUCMD, *) NREPS|\newline
\verb|READ (LUCMD, *) (IDOREP(I), I = 1, NREPS)|

Consider perturbations of selected
symmetries only.  Read one more line specifying how many
symmetries, then one line listing the desired symmetries. Note that
only those symmetries previously defined to be true with the keyword
\Key{REPS} from the \aba\ input modules will be calculated. This
keyword thus represents a subset of the \Key{REPS} of the general
input module.

\item[\Key{RESTART}] Tells the program that this is a restarted
geometry optimization\index{restart!geometry optimization} and that
information may therefore be available
on the \verb|DALTON.WLK| file.

\item[\Key{REUSE}] Use the property derivatives available on the file
\verb|DALTON.WLK| in a calculation of the harmonic contribution to the
vibrational average. In this case, only a new force field will be calculated.

\item[\Key{SCALE}]\verb| |\newline
\verb| READ (LUCMD, *) NUMNUC|\newline
\verb| DO INUC = 1, NUMNUC|\newline
\verb|    READ (LUCMD, *) IATOM,(SCALCO(J,IATOM), J = 1, 3)|\newline
\verb| END DO|

Scale the atomic coordinates.  Read one more
line specifying how many atoms to scale, then one line for
each of these atoms~(*) specifying the atom number and scale
factors for all three Cartesian coordinates. Default is no scaling of
the atomic coordinates.

%\item[\Key{STRICT}] Strict mode following. Obsolete keyword. Do not use.

\item[\Key{TEMPERATURES}]\verb| |\newline
\verb|  READ (LUCMD, *) NTEMP|\newline
\verb|  READ (LUCMD, *) (TEMP(ITMP), ITMP=1,NTEMP)|

Read a set of temperatures for which the
effective (rovibrationally averaged) geometries are to be
calculated. Read one more line containing the number of different
temperatures, and another line containing the list of temperatures.

\item[\Key{TOLERANCE}]\verb| |\newline
\verb|READ (LUCMD, *) TOLST|

Threshold for convergence of the geometry optimization (on gradient
norm).  Read one more line specifying the threshold~(*).  Default
is~10$^{-5}$.

\item[\Key{TRUST}]\verb| |\newline
\verb|READ (LUCMD, *) TRUSTR, TRUSTI, TRUSTD|

Trust region information\index{trust radius!.WALK geometry optimization}.  Read one more
line specifying three values~(*): initial trust radius, factor by
which radius can be incremented, and factor by which it can be
decremented.  Defaults are~0.5, 1.2 and 0.7, respectively; initial
trust radius default is~0.3 if desired molecular Hessian index is
greater than zero. In dynamic walks\index{dynamics} the trust radius
is by default put
to 0.005, and in walks along an Intrinsic Reaction Coordinate
(IRC)\index{intrinsic reaction coordinate}\index{IRC} the
default trust radius is 0.02. For dynamical walks the default
increment and decrement factor is changed to~2.0 and~0.8
respectively.

\item[\Key{VIBAVE}]
Request the calculation of the harmonic contribution to the
vibrational average of a molecular property.

%\item[\Key{VIBCNV}] Requests a vibrational analysis at the optimized
%geometry. Obsolete keyword. Do not use.

\item[\Key{ZERGRD}]\verb| |\newline
\verb|READ (LUCMD, *) ZERGRD|

Threshold below which gradient elements are
treated as zero.  Read one more line specifying value~(*). Default
is~10$^{-5}$. This keyword is mainly used for judging which modes are
symmetry breaking when using the keyword \Key{KEEPSY} as well as
when deciding what step to take when starting a walk from a transition
state.
\end{description}

\subsection{Molecule geometry and basis sets, \Sec{MOLBAS}}\label{sec:herrdn}

The directives in this section extend or modify the reading of the molecule geometry and basis set
specifications in the "\verb|.mol|" input file (internally \molinp).

\begin{description}
\item[\Key{CM FUN}]\verb| |\newline
\verb|READ (LUCMD,*) LCMMAX, CMSTR, CMEND|

Use Rydberg basis functions\index{Rydberg basis function}\index{basis function!Rydberg}
(center of mass functions\index{center of mass function}) as suggested by
Kaufman {\it et al.\/}~\cite{kkwbmjjpbamop22}. \verb|LCMMAX| denoted
the maximum quantum number of the Rydberg functions, basis functions
for all quantum up to and including \verb|LCMMAX| will be generated
(s=0, p=1 etc.) \verb|CMSTR| and \verb|CMEND| are the half-integer
start- and ending quantum number for the Rydberg basis functions. The
basis functions will be placed at the position of a dummy center
indicated as \verb|X| in the \molinp\ file, for example
\begin{verbatim}
Charge=0.0 Atoms=1 Basis=CMFUN
 X    0.0      0.0       0.0
\end{verbatim}
The charge of
the ion-core is determined by the keyword \Key{ZCMVAL} (default: +1),
not by a charge specified for \verb|X|; use \verb|Charge=0.0|.
If no center named \verb|X| is present in the \molinp\ file,
this input will be ignored.
(If you use \verb|ATOMBASIS| in the \molinp\ file, you still need to specify
a basis set name for the center \verb|X|, as in the example. However, you can write anything
because it will be ignored. Maybe for clarity you would like to use \verb|Basis=CMFUN|.)


\item[\Key{CONTINUUM}]\verb| |\newline
\verb|READ (LUCMD,*) LCNTMAX, CNTSTR, CNTEND|

%
% Introduced by Sonia Coriani, 2012
%
Use continuum-like basis functions\index{Continuum basis function}\index{basis function!Continuum}
(Gaussian continuum-like\index{Gaussian continuum basis}) as suggested by
Kaufman {\it et al.\/}~\cite{kkwbmjjpbamop22} (cf. eq. 20). \verb|LCNTMAX| denoted
the maximum quantum number of the continuum functions, basis functions
for all quantum up to and including \verb|LCNTMAX| will be generated
(s=0, p=1 etc.) \verb|CNTSTR| and \verb|CNTEND| are the integer
start- and ending quantum number for the continuum basis functions. The
basis functions will be placed at the position of a dummy center
indicated as \verb|X| in the \molinp\ file.
See example above for the \Key{CM FUN} input.
No charge of the ion-core is required for these functions.
If no center named \verb|X| is present in the \molinp\ file,
this input will be ignored.

\item[\Key{MAXPRI}]\verb| |\newline
\verb|READ (LUCMD,*) MAXPRI|

Set maximum number of primitives\index{primitive orbitals!maximum number} in any
contraction.  Read one more line containing number.  Default
is~25, except for the Cray-T3D/E, where the default is 14.

%\item[\Key{MOLINP}] Indicates that the molecular input comes at the
%end of the current file.
%Default is that the molecular input is to be read from unit~9.

\item[\Key{NOMOVE}]\verb| |\newline
Do not move molecule to center-of-mass if doing automatic symmetry detection.
Do not rotate the molecule after automatic symmetry detection to get maximum symmetry in the calculation
(a rotation axis or a mirror plane normal vector must be along a coordinate axis for \dalton\ to use it).

\item[\Key{NUCMOD}]\verb| |\newline
\verb|READ (LUCMD, *) INUC|

Choose nuclear model.
A 1 corresponds to a point nucleus (which is the default in Dalton),
and 2 corresponds to a Gaussian distribution model (which is the default in the Dirac program, \url{http://diracprogram.org}).

\item[\Key{PRINT}]\verb| |\newline
\verb|READ (LUCMD,*) IPREAD|

Set print level for input processing.  Read one more line containing
print level. Default is the \verb|IPRDEF| from the \Sec{*INTEGRALS} input
module.

\item[\Key{R12AUX}] An auxiliary basis is used.
Basis sets must be identified as either orbital basis or
auxiliary basis in the \molinp\ file (line 5).
See Sec.~\ref{sec:r12aux} on page~\pageref{sec:r12aux}.

\item[\Key{SYMTHR}]\verb| |\newline
\verb|READ (LUCMD,*) TOLLRN|

Read threshold for considering atoms to be related by symmetry. Used
in the automatic symmetry detection routines. Default is $5.0\cdot
10^{-6}$.

\item[\Key{UNCONT}]

Force the program to use all input basis sets as primitive (completely decontracted) sets.

\item[\Key{WRTLIN}]

Write out the lines read in during the input
processing of the \molinp\ file. Primarily for debugging
purposes or for analyzing input errors in the \molinp\ file.

\item[\Key{ZCMVAL}]\verb| |\newline
\verb|READ (LUCMD,*) ZCMVAL|

Read the charge of the ion-core center for the Rydberg basis functions specified
by the \Key{CM FUN} keyword. Default is a charge of one.

\end{description}

\section{Numerical differentiation : \Sec{*NMDDRV}}\label{sec:nmddrv}

This module can calculate any geometrical derivative of the energy
using either as high analytical derivatives as possible, or using the
specified level of analytical derivatives (assuming implemented for
the choice of wave function)~\cite{numder}. Also performs vibrational averaging over
selected first- and second-order molecular properties.

\begin{description}

\item[\Key{DISPLACEMENT}]\verb| |\newline
\verb|READ (LUCMD,*) DISPLC|

Reads in the step lengths (in atomic units) that is to be used in the
numerical differentiation scheme. Default is 1.0D-2.
This value is also used for any numerical gradients and Hessians for
geometry optimizations controlled by the \Sec{OPTIMIZE} module.

\item[\Key{DORDR}]\verb| |\newline
\verb|READ (LUCMD,*) NMORDR, NAORDR|

Sets the numerical (NMORDR) and analytical (NAORDR) differentiating
order for calculating force constants. Current implementation has an
artificial boundary at the 5. numerical derivative (independent of the
analytical differentiation order). Notice that if you would like to
calculating 4.derivatives from analytical 2.derivatives the input
would be \verb|2  2|, since you would like to get 2.order numerical
derivatives from 2.order analytical derivatives.

\item[\Key{DRYRUN}]\verb| |\newline
\verb|READ (LUCMD,*) NMREDU|\newline
\verb|READ (LUCMD,*) (KDRYRN(II),II=1,NMREDU)|

The numerical derivatives will not be calculated, and the program will
just set up the required displacement. An optional number of redundant
coordinate displacements can be specified, corresponding to translational and
rotational degrees of freedom.

\item[\Key{HARMONIC FORCE FIELD}]\verb| |\newline
At the end of the calculation, perform a vibrational analysis using
the calculated molecular Hessian matrix.

\item[\Key{MANUAL}]\verb| |\newline
Dump each individual geometry input file to the \verb|DALTON.OUT|
file. Primarily of interest for debugging purposes.

\item[\Key{NORMAL}]\verb| |\newline
Using this keyword the numerical differentiation will be carried out
with respect to normal coordinates. The program will then do the
necessary numerical differentiation to get the molecular Hessian (and thus
normal coordinates), before it will carry on calculating the
higher-order force field requested with respect to the calculated normal
coordinates.

\item[\Key{PRECALCULATED HESSIAN}]\verb| |\newline
Use a precalculated molecular Hessian available on the \verb|DALTON.HES| file
when defining normal coordinates. Only active in combination with the
keyword \Key{NORMAL}.

\item[\Key{PRINT}]\verb| |\newline
\verb|READ (LUCMD,*) IDRPRI|

Control the print level in the numerical derivative routines. Default
is the same as the general print level \verb|IPRUSR|.

\item[\Key{PROPAV}]\verb| |\newline
Indicates that an averaging over vibrational motions of a molecular
property (including the energy) is to be performed. The input
for what kind of a vibrational analysis is to be performed is
specified in the \Sec{PROPAV} input submodule.

\item[\Key{PROPER}]\verb| |\newline
\verb|READ (LUCMD,*) NMRDRP, NARDRP|

Sets the numerical (NMRDRP) and analytical (NARDRP) differentiating
order for calculating geometrical derivatives of molecular
properties. Currently, the possible presence of analytical property
derivatives cannot be taken advantage of, and the default value of $0$
for \verb|NARDRP| should be used. Currently, dipole transition
strengths and vibrationally averaged spin--spin coupling constants
have been implemented, the former also for Coupled--Cluster wave
functions\index{CC}\index{Coupled Cluster}. Notice that vibrationally
averaged spin--spin coupling constants also can be calculated using
the \Key{VIBANA} and \Sec{VIBANA} keywords and input section.

\item[\Key{RESTRT}]\ \newline
Controls the restart procedure. If a calculation crashes during a
force constant calculation, there will be a file in the work directory
called RSTRT.FC. This will contain all of the information necessary to
restart the calculation. If this file is available and the keyword is
used, \dalton\ will attempt to restart the calculation.

\item[\Key{REUSE HESSIAN}]\verb| |\newline
If a (mixed) numerical molecular Hessian has been calculated, it will be saved
in the file \verb|DALTON.HES| for possible future use.

\item[\Key{SDRTST}]\verb| |\newline
If analytical molecular Hessian also has been requested and second-order
numerical derivatives, a comparison of the numerical and analytical
molecular Hessians will be done. Primarily used for testing.

\item[\Key{SPECTRO INTERFACE}]\verb| |\newline
Write an interface file \verb|DALTON.SPC| containing force fields of
different orders suited for analysis with the SPECTRO
program~\cite{spectroref}.

\item[\Key{SYMMET}]\verb| |\newline
\verb|READ (LUCMD,*) FCLASS|

Assigns the molecular point group of the molecule (FCLASS). For
instance for water FCLASS would be equal to ``C2v''. Main rotational
axis needs to be set to the z-axis in the .mol file. Additional
generating elements needs to be the x-axis for a C2 rotation, and the
xy plane for a mirror plane. The current implementation only allows
for symmetry use, when differentiating from energies.

\item[\Key{TEST N}]\verb| |\newline
Test if the normal coordinates used for calculated geometrical
derivatives give rise to force fields with appropriate
symmetries. Mainly for debugging purposes.

\end{description}

\subsection{Vibrational averaging of molecular properties: \Sec{PROPAV}}
\label{sec:nmddrv.propan}

This module sets up overall numerical or mixed-numerical/analytical
geometrical property derivative calculations, as well as performs selected
post-analyses of the calculated energy and property derivatives in
terms of zero-point vibrationally averaged properties.

\begin{description}

\item[\Key{ANHA-P}]\verb| |\newline
Requests the calculation of the anharmonic contribution to a
vibrationally averaged molecular property, requiring the use of the
first-order perturbed vibrational wave function (requiring knowledge
about parts of the cubic force field) and the first derivative of the
molecular property. See also the keyword \Key{HARM-P}. Currently only
implemented for spin--spin coupling constants.

\item[\Key{CUBIC }]\verb| |\newline
Calculate the molecular Hessian and third derivatives and write the files
\verb|DALTON.HES| and \verb|DALTON.WLK| for use in future
calculation of the vibrationally averaged properties using the
\Key{ONLY-P} keyword.

\item[\Key{EFFECTIVE GEOMETRY}]\verb| |\newline
The effective geometry is calculated. This geometry corresponds to the
zero-point vibrationally averaged geometry and is thus often referred
to as the $r_z$ geometry. For more information on the effective
geometry, see Refs.~\cite{poakrprtjcp112,krpoaprtjcp112}.

\item[\Key{HARM-P}]\verb| |\newline
Requests the calculation of the harmonic contribution to a
vibrationally averaged molecular property, requiring the use of the
unperturbed vibrational wave function and the second derivative of the
molecular property. See also the keyword \Key{ANHA-P}. Currently only
implemented for spin--spin coupling constants.

\item[\Key{MODE ANALYSIS}]\verb| |\newline
Request an analysis of the contribution of the different vibrational modes
to the total zero-point vibrational corrections to a molecular property.

\item[\Key{ONLY-P}]\verb| |\newline
Calculate only the equilibrium value and derivatives of the molecular property
required to calculate the (harmonic and anharmonic) vibrational average.
In order for this to work, \verb|DALTON.HES| and \verb|DALTON.WLK| files containing the molecular Hessian
and third derivatives of the molecule must be available from a previous calculation.

\item[\Key{P-BASIS}]\verb| |\newline
\verb|READ (LUCMD,*) PRPBTX|

When calculating vibrationally averaged properties use the basis set in the
.mol file for calculating the force-constants and \verb|PRPBTX| for
calculating the property and property derivatives.

\item[\Key{SPIN-SPIN COUPLINGS}]\verb| |\newline
Requests the calculation of indirect spin--spin coupling constants.
If combined with the keywords \Key{HARM-P} and \Key{ANHA-P}, the
zero-point vibrational corrections to these coupling constants will be
calculated.

\end{description}

\subsection{Vibrational analysis: \Sec{VIBANA}}
\label{sec:nmddrv.vibana}

This section is identical to the vibrational analysis module described
in Chapter~\ref{sec:abavib}, but appears in the \Sec{*NMDDRV} module
if the molecular Hessian has been calculated numerically using the \Key{NMDDRV}
keyword. We refer to Chapter~\ref{sec:abavib} for a list of available
keywords.
